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Abstract 

In the present study, we identify the need for high accuracy 
computing for capturing the transition to turbulence for the 
flow past a natural laminar flow aerofoil by direct numerical 
simulation. The major computational elements which 
contribute error in performing DNS have been identified 
here and systematically shown to be removed by 
appropriate grid, numerical method and post-processing 
tools. We have shown a better way of generating orthogonal 
grid around thick and cambered aerofoil which helps in 
significant reduction of numerical errors. DNS of SHM-1 
aerofoil for Re = 10.3 * 10 6 has been made possible by the use 
of high accuracy dispersion relation preserving compact 
scheme, which also has the potential to introduce aliasing 
error and spurious upstream propagating disturbances (q- 
waves). Systematic analysis of such methods helps us to use 
an appropriate fine grid along with upwind and multi- 
dimensional filters to remove these problems. The 
simulation results clearly identify the route of bypass 
transition for this flow field. 
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Introduction 

Performance of an aircraft can be improved by a well- 
designed wing. In this context, flow over an aerofoil is 
an important research topic in aerodynamics. The 
present work is focused on obtaining time accurate 
numerical solutions for transitional flows past a 
natural laminar flow (NLF) aerofoil. By definition, 
NLF aerofoils are designed in such a way that 
transition from laminar to turbulent flow is delayed on 
the suction surface, as far aft as possible. The process 
of transition begins via flow instability and in most of 
the design efforts of NLF aerofoils, the onset of flow 
instability is delayed. One of the most commonly 
studied and understood flow instability for wall 
bounded flows is the classical route associated with 
the formation of Tollmien-Schlichting (TS) waves. 



Early researches in flow instability studies focused 
attention on understanding basic mechanisms for 
inviscid disturbance flows, with the assumption that 
viscous actions can only stabilize flows. This is the 
basis of early works of Kelvin, Helmholtz, Rayleigh 
and others, which can be found in (Drazin and Reid, 
1981; Sengupta, 2012). This could not, however, explain 
the flow transition over a flat plate and led to efforts 
which attempted to understand the role of viscous 
actions via the solution of Orr-Sommerfeld equation, 
through the researches reported in (Heisenberg, 1924; 
Schlichting, 1933; Tollmien, 1931), where the stability of 
flow over a flat plate was investigated as an example 
of zero pressure gradient Blasius boundary layer. It is 
shown that if small amplitude monochromatic time 
harmonic disturbance is excited locally, then the 
response field appears as TS wave and travels in space 
exhibiting growth. Results in (Schlichting, 1933; 
Tollmien, 1931) predicted onset of criticality and 
provided properties of unstable waves, whose 
experimental verification came via the pioneering 
effort of (Schubauer and Skramstad, 1947). Since then 
spatial growth of TS waves is thought to be 
responsible for flow transition with the implicit view 
that the primary instability eventually leads to flow 
transition, following other secondary and tertiary 
instabilities. While different frequencies suffer 
differential growth rates, appearance of TS wave is 
considered as a precursor to transition. 

According to (Morkovin, 1969), there are instances 
where flow transitions are seen to occur without the 
presence of TS waves and these are collectively 
categorized as bypass transition. We note that such 
terminology does not imply that there is only one 
bypass transition route. For example, one such bypass 
transition route has been shown theoretically and 
experimentally in (Lim et ah, 2004) and (Sengupta et ah, 
2003). The physical mechanism has been identified 
with the help of an energy based receptivity equation 
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obtained from Navier-Stokes equations without 
making any approximation. An example of this route 
of bypass transition on the attachment line, along the 
leading edge of an infinite swept wing, has been 
shown via convecting aperiodic vortices in the 
freestream as sub-critical instability (Sengupta and 
Dipankar, 2005). 

The classical or the bypass route of transition has been 
explained above for canonical flows only, while it is 
known that adverse pressure gradient has strong 
destabilizing effects on flow transition. As NLF 
aerofoils for general aviation purposes are of large 
thickness, the aft portion of the aerofoil experiences 
very severe varying pressure gradients, on both the 
top and the bottom surfaces for different angles of 
attack. Existence of such varying adverse pressure 
gradient is capable of amplifying very small 
background disturbances and no systematic criterion 
exists which theoretically correlates background 
disturbances with transition to turbulence for either 
the classical or the bypass route. In this context, it is 
worthwhile to investigate the process of transition 
over an aerofoil by direct numerical simulation of 
Navier-Stokes equation. 

Here, we have studied flow past a NLF SHM-1 
aerofoil (Fujino et al, 2003) which has been designed to 
delay flow transition at a cruise Reynolds number of 
20.3 million. This is usually done by contouring the 
section, so that the pressure gradient delays flow 
transition on the suction surface, as far as possible. 
However, when such an aerofoil is employed for 
general aviation purpose, the flow eventually suffers 
transition in the rear half of the section due to the 
prevalent adverse pressure gradient which accelerates 
the growth rate of the trace amount of background 
disturbances. In a computational framework, the 
background disturbances can arise from various 
sources of numerical errors. In nature, disturbances 
are omnipresent and their role in destabilizing fluid 
flow is a pacing item of research in fluid dynamics. 
These disturbances can arise from surface 
perturbations, free-stream noise, surface vibrations etc. 
We emphasize that even when explicit excitation is not 
employed in numerical computations, due to 
hypersensitivity of the flow field to background 
disturbances at higher Reynolds numbers, transition is 
observed. In an ideal scenario, one would like to 
remove all sources of numerical errors which 
contribute to flow transition and then apply explicit 
excitations which mimic the physical background 
disturbances. This is possible for Blasius boundary 



layer and results are reported in (Fasel and Konzelmann, 
1990; Sengupta and Bhaumik, 2011; Sengupta et al, 2009a). 
It is therefore imperative to understand the various 
sources of numerical errors and the way these can be 
eliminated or reduced. It is also equally important to 
choose an appropriate formulation and its accurate 
discretization. 

The sources of errors while discretizing are (i) aliasing 
error in computing product terms, as explained in 
(Sengupta, 2004; Sengupta et al, 2009b); (ii) spurious 
dispersive waves (whose extreme form is known as q- 
waves), as described in (Poinsot and Veynante, 2005; 
Sengupta et al, 2012a; Vichnevetsky, 1982) and due to 
stability, phase and dispersion error, as explained in 
(Sengupta et al, 2007). 

This paper is formatted in the following manner. In 
the next section, we provide the governing equation 
for 2D DNS performed here and its formulation. In 
section 3, we describe an improved orthogonal grid 
generation method which helps in significantly 
reducing numerical error. Numerical methods and 
analyses of these are provided in section 4. Detailed 
results of DNS past the SHM-1 airfoil are provided in 
section 5. The paper closes with a summary and 
conclusions. 

Governing Equations and Formulation 

We have obtained numerical solution of Navier-Stokes 
equation using stream function (ip) - vorticity (co) 
formulation. This formulation ensures solenoidality of 
velocity and vorticity field simultaneously and found 
to be more accurate, as compared to other 
formulations (Sengupta, 2004). Also, we employ an 
orthogonal grid for the studied incompressible flow 
which allows explicit use of orthogonal formulation 
resulting in fewer numerical computations per time 
step and thereby reducing numerical errors by 
discretizing fewer terms with efficient discretization 
following the physical nature of the flow. This latter 
aspect is best demonstrated by the isotropic treatment 
of the diffusion operators in orthogonal formulation, 
as compared to non-orthogonal formulations. The 
governing Navier-Stokes equation is decoupled into 
kinematics and kinetics of the flow in (ip, co)- 
formulation and expressed in an orthogonal grid by 
the following equations 
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where hn and h22 are the scale factors used in mapping 
physical (x, y)-plane to a computational (£, 7])-plane 
(Sengupta, 2004), where E, co-ordinate is in azimuthal 
direction and r\ is normal to it. For orthogonal 
mapping (Nair and Sengupta, 1998; Sengupta et al, 2007), 
one can define the scale factors by, 



and 



The leading and trailing edge portions of the aerofoil 
have a sharp slope variation and to simulate the flow 
correctly distribution of points on the input aerofoil 
geometry has to be performed judiciously by 
clustering more number of points in those regions. In 
order to resolve the boundary layer formed over the 
aerofoil, we have clustered the points in the r\ direction 
using a tangent hyperbolic distribution given later. 

Following initial and boundary conditions are used in 
solving the governing equations. On the aerofoil 
surface, following no-slip conditions are used 



.(3) 



iff = constant, = 

dn 



These conditions also help fix the wall vorticity, which 
is required as the boundary condition for the vorticity 
transport equation, Eq. (2). In O-grid topology, one 
introduces a cut starting from the leading edge of the 
aerofoil to the outer boundary. Periodic boundary 
conditions apply at the cut, which are introduced to 
make the computational domain simply-connected. 
For the stream function equation, Eq. (1), at the outer 
boundary, Sommerfeld boundary condition is used on 
the ^-component of the velocity field. Resultant value 
of streamfunction is used to calculate the vorticity 
value at the outer boundary. 

From the stream function equation the wall vorticity is 
calculated as 



..(4) 



- body 



An Improved Orthogonal Grid Generation 
Method 

For an orthogonal mapping from physical to 
transformed plane, the included angle between the £, = 



constant and 77 = constant lines in the physical plane is 
constrained by 



•(5) 



where the subscripts indicate partial derivative with 
respect to these. One can introduce a distortion 
function / given by 

/ /• (6) 

Using Eq. (5) in Eq. (6), one gets the following two 
relations 

y n =-fi [ i ( ? ) 

*„ = fy e (8) 

These are the well-known Beltrami equations for 
general orthogonal mapping (Duraiswami and 
Prosperetti, 1992; Nair and Sengupta, 1998; Ryskin and 
Leal, 1983; Sengupta, 2004) and these are nothing but 
the generalization of Cauchy-Riemann relations. The 
authors in (Dipankar and Sengupta, 2006; Nair and 
Sengupta, 1998) have used hyperbolic grid generation 
method to maintain orthogonality in a stricter sense 
for generated grids. This has been achieved by using 
Eq. (5) as one of the governing equations. In (Steger and 
Chaussee, 1980), prescription of the cell volume given 
by hnh22 has been used as the second governing 
equation for outward marching of the equations from 
the boundary. In (Nair and Sengupta, 1998), one of the 
Beltrami equations, Eqs. (7) and (8), has been used 
along with the orthogonality condition given by Eq. 
(5). 

In the present work, an improved procedure for grid 
generation around a SHM-1 aerofoil is described step- 
by-step. The SHM-1 aerofoil has a camber line which 
results in concave surface near the trailing edge on the 
bottom surface, as shown in the top frame of Fig. 1(a). 

In this figure, aerofoil is shown with a box A-B-C-D 
marking the concave part of the trailing edge. If the 
procedure given in (Nair and Sengupta, 1998) is 
followed up to a short distance near the aerofoil 
surface, then the grid-lines in the region A-B-C-D are 
as shown in the bottom frame of Fig. 1(a). However, if 
this process is continued further, then the grid-lines 
will eventually intersect, as indicated by the 
converging arrowheads in the figure. We note that 
even when the grid-lines move towards each other 
without intersection, the metric variation causes 
numerical difficulties associated with solution of 
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Navier-Stokes equation. Thus, it is imperative that 
grid-lines should emerge smoothly from the surface 
without the possibility of intersection. 




SHM- 1 aerofoil with a region ABCD marked 
near the trailing edge showing concave surface 



Zoomed view of the generated grid in 
the region ABCD after Step I. 




Terminal line 



FIG. 1(a) SCHEMATIC OF SHM-1 AEROFOIL IS SHOWN IN 
THE TOP FRAME (I) WITH THE MARKED REGION ABCD 
HIGHLIGHTING CONCAVE PORTION OF THE AEROFOIL. 
BOTTOM FRAME (II) SHOWS ZOOMED VIEW OF THE 
GENERATED ORTHOGONAL GRID USING THE PROCEDURE 
GIVEN IN (NAIR AND SENGUPTA, 1998). NOTE THE ARROWS 
ARE POINTED TOWARDS EACH OTHER; INDICATE PROBLEM 

OF GRID INTERSECTION IF THE GRID IS GENERATED 
CONTINUOUSLY FOLLOWING THE METHOD IN (NAIR AND 
SENGUPTA, 1998). 

The problem is acute due to the presence of concave 
region on the basic geometry or in the generated grid- 
lines subsequently. Grid intersection problem 
associated with concave surface can be avoided as 
explained below. 

(I) Create an orthogonal grid around a SHM-1 aerofoil 
for a limited region close to the surface as shown in 
bottom frame of Fig. 1(a) following the procedure in 
(Nair and Sengupta, 1998). 

This is reported here as Step I and the generated grid 
consists of 597 points in the azimuthal direction and 80 
points in the t] direction. The grid is generated by 
taking identical spacing on r\ = constant line, up to 
0.025c distance from the surface of the aerofoil, where 
c is the chord of the aerofoil. The gap between 
successive r\ = constant lines is stretched in the r\ - 
direction following the distribution given by 



K 2 =H 



( tanh [/?(!- 77)] 
tanh[/?] 



.(9) 



where /3 represents the clustering parameter. The grid 
spacing in r/-direction is represented by grid scale 
parameter hn and H is the non-dimensional (in terms 
of chord of the aerofoil) distance of the outer boundary. 



(II) To circumvent problems of grid intersection, 
concavity is gradually removed from successive r/ = 
constant lines as one marches outward from the 
aerofoil. The concave region is identified by plotting 
d 2 y/dx 2 parameter as a function of points on the 
terminal line of the grid region. The concave part of 
the terminal line is shown clearly by the positive 
values on the top surface and negative values below 
the bottom surface. For the generated grids' terminal 
line, this is shown dotted in Fig. 1(b). 

If I122 in the concave region is increased, then the 
concavity will be reduced progressively as we move 
away from aerofoil surface, as shown by the circular 
symbols in Fig. 1(b). We have limited the absolute 
magnitude of the curvature on top and bottom 
surfaces to a threshold value of 0.25 to ensure removal 
of concavity. Threshold value is reached by changing 
hn by a small magnitude in an iterative fashion, as 
rapid change causes grid intersection. Although, the 
problem associated with grid line intersection due to 



Curvature of a terminal line over 
the top surface of aerofoil. 
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Curvature after removal of concavity 



0.2 



0.4 



x/c 



0.6 



0.8 



(in 



X 
33 



o 
o 
o 
o 
o 
o 
o 
o 
o 
o 
o 
o 
o 
o 



Curvature of a terminal line below 
the bottom surface of aerofoil. 



Curvature after Step I 

Curvature after removal of concavity 




Negative values 
show 
concavity 



0.5 

x/c 



FIG. 1(b) VARIATIONS OF CURVATURE ON THE TERMINAL 
LINE OF THE GENERATED GRID ARE SHOWN: (I) OVER THE 
TOP SURFACE OF THE AEROFOIL AND (II) BELOW THE 
BOTTOM SURFACE OF THE AEROFOIL. NOTE THAT THE 
CONCAVITY FROM THE TERMINAL LINE OF THE GENERATED 
GRID HAS BEEN REMOVED COMPLETELY. 
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concavity is removed, one ends up with curvature 
discontinuities, as marked by points At, Ai and Bi, B2 
in Fig. 1(b). Such discontinuities cause discontinuous 
variation of grid metrics. This is avoided by the 
following. 

Next, we do not prescribe a constant threshold limit at 
all points on the terminal line, instead the threshold 
value varies with These threshold limits are 
obtained by using weighted average procedure at the 
i th node by 

d 2 y ( d 2 y d 2 y ) 

—irl =0.50 ^ L, + !L t\, , 
dx 2 ' [dx 2 1+1 dx 2 "'J 

The h.22 is again adjusted in an iterative manner so that 
points on the terminal line achieve the prescribed 
smooth curvature. This procedure removes the 
discontinuities in the curvature, as shown in Fig. 1(c). 
As noted earlier, curvature was discontinuous as 
marked by points Ai, Ai and Bi, Bi. This step removes 
the sharp curvature variation with addition of few 
more points, as shown in Fig. 1(c). 

(i) Curvature of a terminal line over 

the top surface of aerofoil. 

o Curvature after removal of concavity 
, Curvature_after additonal_smc«fJiing 




3 — 1 — 1 — 1 — ois — 1 — 1 — 1 — r 

x/c 



FIG. 1 (c) VARIATION OF CURVATURE ON THE TERMINAL 
LINE OF THE GENERATED GRID IS SHOWN: (I) OVER THE TOP 
SURFACE OF THE AEROFOIL AND (II) BELOW THE BOTTOM 
SURFACE OF THE AEROFOIL AT THE END OF STEP E 



Grid generated at the end of this step is shown in the 
top frame of Fig. 1(d). Marked arrows clearly show 
that the grid concavity has been removed, as also seen 
in the bottom frame of Fig. 1(d). Here, we have shown 
the trailing edge of SHM-1 aerofoil along with the 
terminal lines of the generated grid after Steps I and II. 
Concavity is thus, removed completely by the 
procedure of filling the gap in Step II. 

(Ill) In this step, we further generate grid lines by 
marching away from the terminal line of Step II up to 
desired outer boundary. This step completes the grid 
generation process. 




FIG. 1(d) ZOOMED VIEW OF THE GENERATED 
ORTHOGONAL GRID AT THE END OF STEP II IS SHOWN IN 
THE TOP FRAME (I). NOTE THE CONCAVITY PRESENT IN THE 
INITIALLY GENERATED GRID HAS BEEN REMOVED AT THE 
END OF STEP II, AS SHOWN IN THE BOTTOM FRAME (II) 



Zoomed view of the generated grid around SHM- 1 aerofoil 




(b) 
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deviation from orthogonality A9 
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FIG. 2(a) ZOOMED VIEW OF THE GENERATED GRID AFTER 
STEP III IS SHOWN; (b) CONTOURS OF DEVIATION FROM 
ORTHOGONALITY ARE SHOWN FOR THE GENERATED GRID 
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Here we have constructed a final grid of 597 X 397 
points and outer boundary located at a distance of 14c 
from the aerofoil surface. A zoomed view of the 
generated grid is shown in Fig. 2(a). Deviation of grid 
lines from orthogonality is shown in the contour plot 
of Fig. 2(b). The maximum deviation from 
orthogonality {AOmax) is 0.0748° and the minimum 
deviation (Admin) is -0.0424°. This is one order less 
deviation as compared to 0.524° reported for the case 
of NACA 0015 aerofoil in (Nair and Sengupta, 1998). 
This establishes the present grid generating algorithm 
that has a capability to generate grid around concave 
shaped bodies with negligible deviation from 
orthogonality. 

Analyses of Numerical Methods 

Before obtaining numerical solutions, we show 
numerical properties of the used high accuracy space- 
time discretization scheme (SOUCS3-OCRK3), as 
explained with properties in Fig. 3. Details of the space 
and time discretization methods can be seen in 
(Sengupta et ah, 2011a). Firstly, we compare the spectral 
resolution of SOUCS3 method with some conventional 
spatial discretization schemes used to obtain the first 
derivatives in Eq. (2). 

If a variable is defined by its Laplace transform as 
u(x,t) =J U(k,t) e ikx dk, then its exact spatial derivative is 



■ jik U(k,t)e' b dk 



When this derivative is obtained numerically, it is 
written equivalently as ' ^ u 



-\ik ei U(k,t)e ik 



dk 



In Fig. 3(a), spatial resolution for the first derivative 
given by real part of keq/k for central second order 
scheme CDr, for fourth order central discretization 
scheme CDi, along with the compact scheme SOUCS3 
(Sengupta et al, 2011a), is compared as a function of 
non-dimensional wavenumber (kh), up to the Nyquist 
limit, with h representing the uniform grid spacing in 
the transformed plane. Ideally the ordinate should be 
equal to one across the complete kh range. However, 
phase mismatch amounts to filtering the unknown in 
obtaining the derivative numerically. Note that CDi 
method produces highest filtering at all wavenumbers, 
while SOUCS3 scheme has near-spectral accuracy over 
a larger range. From Fig. 3(a), SOUCS3 is seen to 
resolve first derivative over a range by order of 
magnitude, as compared to the CDi scheme in each 



direction, leading to significantly enhanced ability of 
DNS by such compact schemes. 



SOUCS3 




IGI - contours (SOUCS3+OCRK3) 
(b) Max : 43. 1405, Min : 0.9998 




1 N 2 



V 

4 



/C - contours (SOUCS3 + OCRK3) 
Max : 1.424, Min : -7.54 




1 N 2 



FIG. 3(a) COMPARISON OF SCALE RESOLUTION BETWEEN 
SECOND ORDER, FOURTH ORDER CENTRAL DIFFERENCE 
SCHEME AND THE COMPACT SCHEME SOUCS3 USED HERE; 
(b) I G I AND (c) VgN/C CONTOURS FOR THE SOLUTION OF ID 
CONVECTION EQUATION USING THE SOUCS3-OCRK3 
SCHEME. NOTE THE PRESENCE OF q-WAVE REGION AT 
HIGHER WAVENUMBER RANGE IN (c) 

However, one must also ensure that the adopted 
method correctly propagates energy and phase of 
created disturbance field maintaining perfect neutral 
stability, as emphasized in (Sengupta et al, 2007), with 
respect to the analysis of space-time discretization 
schemes for the linear convection equation 



8t dx 



(10) 



Equation (10) convects the initial solution to the right 
at the group velocity, which is equal to the phase 
speed C, at all time due to non-dispersive nature. 
Following the analysis of (Sengupta et al, 2007), we 
show the numerical amplification factor I G I -contours 
and numerical group velocity VsN/C-contours for the 
solution of Eq. (10) using SOUCS3 and OCRK3 
schemes in Figs. 3(b) and 3(c). 

The numerical stability is defined by the amplification 
factor G (k, At) = U (k, t+At)/U (k, t) which is a complex 
quantity. Numerical dispersion relation is obtained in 
terms of fan/3 = -Gimag/Greai and numerical group velocity 
has been obtained in (Sengupta et al, 2007; Sengupta et 
al, 2011a) as, 

v s n _ i dp 



C N c d(kh) 
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with the Courant-Friedrich Lewy (CFL) number 
defined by Nc = CAt/Ax. 

In Fig. 3(b), IGl-contours are shown for SOUCS3- 
OCRK3 scheme in (Nc, kh)-p\ane for the solution of Eq. 
(10). The IGl-contours have the desirable neutrally 
stable property across the fully resolved length scales, 
only for a narrow range of small Nc values, shown by 
the hatched region in the figure. The other important 
dispersion property of the numerical method relates to 
preserving the actual physical dispersion relation is 
shown in the (Nc, kh)-p\ane, in Fig. 3(c). DRP property 
is retained effectively for a significant kh and Nc range. 
Figure 3(c) also shows that choice of kh should be 
restricted below kh = 2.4, as the solution components 
above this non-dimensional wavenumber travel in the 
opposite to the physical direction, due to negative 
values of numerical group velocities, shown by 
vertical hatched lines in Fig. 3(c). This is the region, 
which will be responsible for spawning q-waves. 
Extreme form of spurious dispersion associated with 
q-waves leads to unphysical propagation which may 
lead to numerical instabilities. 

Use of conventional central one-dimensional (ID) 
filters, as proposed in (Gaitonde et al, 1999) cannot 
rectify the problem of ^-waves for any value of Nc as 
shown in Fig. 3(c). So a new fifth order upwind filter 
has been introduced in (Sengupta et al, 2009), which 
has the additional capability of modifying numerical 
dispersion relation so that q-wave region is removed 
for a small range of Nc. Effectiveness of a filter while 
filtering different wavenumber components is shown 
by its transfer function (Bhumkar and Sengupta, 2011; 
Sengupta et al, 2009) and in Figs. 4(a) and 4(b), we have 
shown real and imaginary parts of the transfer 
functions for the fifth order upwinded filter 
corresponding to a central node. In these plots, we 
have chosen value of filtering coefficient ai = 0.45 and 
upwind parameter r/j = 0.001. Note that transfer 
function of a filter has both real and imaginary parts. 
The real part shown in Fig. 4(a), indicates the amount 
of attenuation imposed by the filter at different scales. 
Simultaneous presence of real and imaginary parts for 
the transfer function of upwind filter ensures 
alteration of the dispersion property. Combined action 
of both real and imaginary parts of the transfer 
function modifies the numerical amplification factor as 
shown in Fig. 4(c) for the solution of Eq. (10), by 
SOUCS3-OCRK3 schemes. 

One observes attenuation of high wavenumber 
components which are responsible for numerical 



instabilities. Corresponding numerical group velocity 
contours are shown in Fig. 4(d) with the unique 
feature of complete removal of q-waves for a small 
range of Nc values up to Nc = Nc- as marked in Fig. 4(d). 
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FIG. 4 VARIATIONS OF (a) REAL PART AND (b) IMAGINARY 
PART OF THE TRANSFER FUNCTIONS OF THE FIFTH ORDER 
UPWIND FILTER (SENGUPTA ET AL., 2009) FOR FILTER 
COEFFICIENT^, = 0.45 AND//, = 0.002 ARE SHOWN ALONG 
WITH MODIFICATIONS IN (c) I G I AND (d) V s n/C CONTOURS 
FOR THE CASE OF FIG. 3. NOTE THE 9- WAVE REGION HAS 
BEEN REMOVED OVER A SMALL RANGE OF N„ IN (d) 

Although, problems associated with ^-waves are 
successfully handled by using ID upwind filter 
(Sengupta et al., 2009), aliasing error is also another 
significant source of error in the computations 
obtained by high accuracy numerical methods. To 
control the aliasing error effectively and at the same 
time to limit the extra computational cost of filtering, a 
multi-dimensional adaptive filter has been introduced 
in (Bhumkar and Sengupta, 2011), where a 2D filter is 
used in a limited region wherever needed. General 
form of 2D filters (Bhumkar and Sengupta, 2011) is given 
by 

u. .+a n ,(u. , ■ . +k , +u ... ) = 

>.J 2/ ^ i-l, j i+l, j i,j-l i,j+l j 

M a 

Z = (" -•" ) ( n ) 

n=0 L 
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where M is the order of the filter and is equal to one 
for a 2D second order filter. The coefficients a„'s are 
related to user defined filtering coefficient aif and their 
role is to smooth selectively the higher wavenumber 
fluctuations. In the above, increase of M by one, 
increases the order of the filter by two. Equation (11) is 
used for points inside the computational domain and 
solved iteratively using Bi-CGSTAB method (Vorst, 
1992). As these equations are solved iteratively, 
significant computational efforts are needed if filtering 
is performed over the complete domain. It has been 
shown in (Bhumkar and Sengupta, 2011) that filtering 
can be performed in an adaptive manner over small 
patches, whenever required. 




a, r = 0.20 




FIG. 5 (a) SCHEMATIC OF THE ADAPTIVE 2D FILTERS 
(BHUMKAR AND SENGUPTA, 2011) WITH THE FILTER SUB- 
DOMAIN SHOWN IN THE TRANSFORMED PLANE. FOR THE 
2D FILTERS, THE TRANSFER FUNCTIONS ARE SHOWN IN 
FRAMES (b) TO (e) FOR DIFFERENT VALUES OF THE FILTER 
COEFFICIENT^. 

For this purpose, one identifies a square sub-domain 
in the transformed plane where filtering is performed 
smoothly, as shown in Fig. 5(a) by the (9 X 9) points 
filter-tent used in the adaptive filter. Variation of filter 
coefficient aif over the tent is also shown in this frame. 
Corresponding transfer function variations in {kin, 
fc;ft/)-plane are shown in frames 5(b) to 5(e). The region 
corresponding to (hh + kjhj) > n contributes to aliasing 



error and the transfer functions of this multi- 
dimensional filter effectively attenuate the aliasing 
error without causing alteration of the lower 
wavenumbers as discussed in (Bhumkar and Sengupta, 
2011). Aliasing error results from the evaluation of 
product terms in a finite grid. Computations of 
nonlinear product terms or even the computations of 
linear terms in a transformed domain lead to aliasing 
error. 

For example in Eqs. (1) and (2), nonlinear product 
terms appear when the computations are performed in 
a transformed plane even for linear diffusion terms. 
This problem accentuates at low Reynolds number in 
Eq. (2) for the viscous diffusion terms, as this is 
multiplied by 1/Re and involves evaluating triple 
product terms. It is to be noted that the product 
evaluations involving derivatives using higher 
accuracy methods are more prone to aliasing error due 
to the inherent filtering operation of derivative by any 
discrete methods excepting Fourier spectrum method. 
However, lower accuracy methods filter more in 
obtaining the derivatives and hence suffer lesser 
aliasing. 



Experimental results 
Numerical results 



Grid size: 597x397 



Angle of attack (a) 1 



(b) 



Re = 10.3x10" 

a = 0° 



C l (Expi.) = 0.2248 Grid size: 597x397 

C, ,„..„.., = 0.1990 
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FIG. 6 (a) NUMERICALLY OBTAINED Cl WITH ANGLES OF 
ATTACK IS SHOWN ALONG WITH EXPERIMENTAL RESULTS 
OF (FUJINO ETAL., 2003); VARIATION OF INSTANTANEOUS O 
AND C,i IS SHOWN IN FRAMES (b) AND (c), RESPECTIVELY; (d) 

VARIATION OF RMS VALUE OF THE AZIMUTHAL 
COMPONENT OF THE VELOCITY OVER AEROFOIL FOR RE = 
10.3 X 10 6 ON A LINE CLOSE TO THE AEROFOIL 



46 



Frontiers in Aerospace Engineering (FAE) Volume 2 Issue \, February 2013 



www.fae-journal.org 



Results and Discussion 

Next, we validate our numerical results with the 
experimental results provided in (Fujino et al, 2003). 
Here, we have reported results for two Reynolds 
number cases, Re = 2.8 X 10 6 and Re = 10.3 X 10 6 
using a (597 X 397) grid. In Fig. 6(a), we have 
compared variation of numerically obtained mean lift 
coefficient (Q) with angle of attack for Re = 2.8 X 10 6 
with the experimental results provided in (Fujino et al, 
2003). Comparison shows a good match between the 
numerical and experimental values. The present 
calculations are strictly two-dimensional for Re = 2.8 X 
10 6 without any explicit forcing for the case shown in 
Fig. 6(a), while the experimental results involve three- 
dimensionality and background tunnel noise. 

Time variations of lift and drag coefficient for 
Re = 10.3 X 10 6 and zero angle of attack are shown in 
frames (b) and (c), respectively. Time averaged values 
of lift and drag coefficients are shown with the 
horizontal dashed lines in these frames. Results show 
a good match between the experimental and 
numerically obtained values. Although, computations 
reported with a (597 X 397) grid show a good 
comparison with the experimental results, next we 
report the numerical results obtained by using a finer 
(5269 X 577) grid specifically for a careful study of 
flow transition over a NLF aerofoil. Variation of RMS 
value of azimuthal component of disturbance velocity 
over top and bottom surfaces is shown in Fig. 6(d) at a 
distance of 0.000056c from the aerofoil surface. 
Variation of RMS component over the aerofoil surface 
shows that after 60% of chord there is sharp rise in 
fluctuations on the top surface indicating flow 
transition similar to the experimental results shown in 
(Fujino et al, 2003). On the bottom surface, one notes 
even a sharper variation of the RMS component of 
azimuthal velocity at x = 0.70c. Correct prediction of 
transition location highlights importance of highly 
space-time accurate numerical solutions obtained here. 

We have computed flow past SHM-1 aerofoil for 
Re = 20.3 X 10 6 at zero angle of attack and Fig. 7(a) 
shows streamfunction contours at t = 4.50 for this flow. 
Top frame shows the flow field around the complete 
aerofoil, while in the bottom frame; flow field near the 
trailing edge of the aerofoil is shown at the same 
instant. Flow near the leading edge of the aerofoil 
experiences acceleration due to the favourable 



pressure gradient. As flow moves towards trailing 
edge, it experiences progressive adverse pressure 
gradient. Due to varying adverse pressure gradient, 
small separation bubbles are formed on the aerofoil 
surface near trailing edge which moves downstream 
along with the flow. These separation bubbles are 
observed in the zoomed view of a trailing edge as 
shown in the bottom frame of Fig. 7(a). As these 
bubbles move towards the trailing edge, these further 
excite the flow field. As shown earlier in Fig. 6(d), flow 
transition on top surface starts after 0.60c, these 
separation bubbles also first appear around the same 
location. In Fig. 7(b), variations of displacement 
thickness b* on top and bottom surfaces of SHM-1 
aerofoil are shown in frame (i) up to 60% of chord of 
SHM-1 aerofoil at t = 4.50. Additionally, variation of a 
steady flow separation parameter used in Falkner- 
Skan analysis, m= (x/Ue )(dUe/ dx), is shown for the top 
and bottom surfaces of SHM-1 aerofoil in frame (ii) at t 
= 4.50. 




(i) 



Top surface 








Bottom surface 






4 











(ii) 





— Top surface 






— Bottom surface 




\ \ 













FIG. 7 (a) STREAM FUNCTION CONTOURS OVER THE 
COMPLETE AEROFOIL AS WELL AS NEAR THE TRAILING 
EDGE ARE SHOWN FOR THE INDICATED PARAMETERS; (b) 
VARIATION OF THE DISPLACEMENT THICKNESS 5* AND 
STEADY FLOW SEPARATION PARAMETER m, ON TOP AND 
BOTTOM SURFACES OF SHM-1 AEROFOIL IS SHOWN IN 
FRAMES (I) AND (II), RESPECTIVELY 
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FIG. 8 (a) VARIATION OF WALL VORTICITY ON THE TOP 
SURFACE OF SHM-1 AEROFOIL HAS BEEN SHOWN FOR THE 
CASE OF RE = 10.3 MILLION AND a = 0°, AT THE INDICATED 
TIMES 



A horizontal dashed line corresponding to steady 
separation flow criterion (m= -0.09) is marked for 
comparison purpose. Figure shows that on the top 
surface near the leading edge, flow experiences higher 
favourable pressure gradient, as compared to the 
bottom surface. Pressure gradient parameter (m) 
smoothly decreases on the top surface after the initial 
peak near the leading edge. However, on the bottom 
surface, pressure gradient parameter remains 
favourable to a larger distance from the leading edge 
of the aerofoil and then it has a sharp variation aft of 
the mid-chord location. 

Variation of surface vorticity on the top and the 
bottom surfaces at different time instants are shown in 
Figs. 8(a) and 8(b), respectively. Due to acceleration of 
flow near the leading edge on the top surface, one 
notices sharp vorticity variation near the leading edge. 
Formation of wavy disturbances on the top surface 
starts approximately from x = 0.55c onwards, while on 
the bottom surface, these disturbances start appearing 
from x = 0.70c onwards. Time variation of surface 
vorticity shows that disturbances on the bottom 



surface have similar range of vorticity variation as 
compared to the disturbances on the top surface. Flow 
field in the second half part of the aerofoil shows 
strong unsteady behaviour, as compared to first half, 
even when no explicit excitation is applied. 
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FIG. 8 (b) VARIATION OF WALL VORTICITY ON THE BOTTOM 
SURFACE OF SHM-1 AEROFOIL HAS BEEN SHOWN FOR THE 
CASE OF RE = 10.3 MILLION AND a = 0°, AT THE INDICATED 
TIMES 



Formation and propagation of separation bubbles on 
the top surface is shown in Fig. 9(a). Variation of the 
azimuthal component of the velocity (u) at a distance 
of 2.242 X 10- 6 c from the top surface of the SHM-1 
aerofoil is traced in this figure, with time. This height 
corresponds to the second azimuthal line (77 = constant 
line) from the aerofoil surface. At f = 0.40, small wavy 
disturbances appear originating from the trailing edge. 
These disturbances move upstream with time as noted 
in the subsequent frames. Presence of adverse 
pressure gradient near the trailing edge magnifies 
these disturbances, which propagate further upstream 
as shown in this figure. The calculations are performed 
on a finer grid with (5269 X 577) points, to resolve the 
length scales correctly and avoid aliasing error. 
Additionally, we have used upwind filter (Sengupta et 
al, 2009) to prevent any (j-waves contaminating the 
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numerical solution. So the upstream propagation of 
the disturbances suggest physical bypass transition 
phenomenon (Sengupta, 2012) and not due to ^-waves. 
Figure 9(a) shows how additional wave-packets are 
introduced upstream of this propagating main 
disturbance, marked as Tv in the frames at t=1.25 and 
later. This induced wave packet (Tv) grows quite 
rapidly, as noted from the frames from t=1.25 to t=1.31, 
while convecting downstream. As time progresses, 
this sequence of upstream migration of disturbance 
fronts fill up the top surface, up to about x = 0.55c, at 
around t =1.50. Thereafter, disturbances formed on the 
top surface stop upstream propagation and continue 
to convect downstream, as shown till t =4.50. 
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FIG. 9 (a) VARIATION OF THE AZIMUTHAL COMPONENT OF 
THE VELOCITY (U) AT A DISTANCE OF 1.242 X10- 6 c FROM THE 
TOP SURFACE OF SHM-1 AEROFOIL HAS BEEN SHOWN FOR 
THE CASE OF RE = 10.3X10<- AND a = 0°, AT THE INDICATED 
TIMES. (CONT) 

Similar formation and propagation of disturbances on 
the bottom surface is shown in Fig. 9(b). In this figure, 
once again the variation of azimuthal component of 
the velocity (u) is shown at a distance of 2.242 X 10~ 6 c 
from the bottom surface. Similar to top surface, wavy 
disturbances originate from the trailing edge and can 
be observed as early as at t = 0.50. By f = 0.67, an 
induced small disturbance packet is noted close to x = 
0.8c, which is marked as B». Its rapid growth is traced 



in the subsequent time frames and one can note the 
nonlinear distortion by t = 0.71. 



' »1 ' ^ u'm ' i; 



FIG. 9 (a) VARIATION OF THE AZIMUTHAL COMPONENT OF 
THE VELOCITY (il) AT A DISTANCE OF 1.242 XlO^c FROM THE 
TOP SURFACE OF SHM-1 AEROFOIL HAS BEEN SHOWN FOR 
THE CASE OF RE = 10.3X10* AND a = 0°, AT THE INDICATED 
TIMES 

The trailing edge position near x = 0.8c of the aerofoil 
has large concavity imposing strong adverse pressure 
gradient, which results in such drastic amplification of 
disturbances. Due to nonlinear action, unsteady 
separation bubbles are formed in this part of the 
aerofoil, resulting in increased drag. As compared to 
the top surface, here the upstream location of the 
disturbance is restricted up to x = 0.70c, as noted in the 
bottom frames of Fig. 9(b). 

Summary and Conclusions 

Here, we report the DNS results of transitional flow 
over a NLF aerofoil (SHM-1) for the cruise Reynolds 
number of 20.3 million. To perform such DNS, one 
must first identify the need for high accuracy 
computing to capture the physical processes of 
transition to turbulence. The major computational 
elements essential for DNS have been identified here 
first by reducing various sources of numerical errors. 
In the present work, we have systematically removed 
these by use of appropriate grid, numerical method 
and post-processing tools. 

An algorithm to generate a truly orthogonal grid for 
thick and cambered aerofoil is developed and 
explained with the help of Figs. 1(a) to 1(d), by 
removing grid shocks arising from the concave part of 
the aerofoil near the trailing edge. The resultant 
orthogonal grid for SHM-1 NLF aerofoil is shown in 
Fig. 2, with highest departure from orthogonality is 
noted for a single point as 0.0748° in Fig. 2(b). Thus, 
one can solve Navier-Stokes equation in orthogonal 
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formulation, which reduces numerical errors as 
compared to non-orthogonal formulation. Efficiency of 
the proposed grid generation algorithm in 
constructing orthogonal grid around cambered 
aerofoils is a significant element of the present 
procedure. 
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FIG. 9 (b) VARIATION OF THE AZIMUTHAL COMPONENT OF 
THE VELOCITY (U) AT A DISTANCE OF 1.242 X 10* c FROM THE 
BOTTOM SURFACE OF SHM-1 AEROFOIL HAS BEEN SHOWN 
FOR THE CASE OF RE = 10.3 X 10 6 AND a = 0°, AT THE 
INDICATED TIMES 



Above mentioned orthogonal grid used with 
orthogonal formulation is further aided by adopting 
high accuracy method in solving vorticity transport 
equation. High resolution ability of the used SOUCS3 
(Dipankar and Sengupta, 2006) scheme in discretizing 
convection terms is demonstrated by comparing the 
equivalent numerical wavenumber in the spectral 
plane with conventional second and fourth order 
discretization methods in Fig. 3(a). The compact 
scheme has higher resolution by almost an order of 
magnitude, as compared to second order 
discretization. This spatial discretization is used to 
develop a DRP scheme for time integration, as the 
OCRK3 scheme in (Sengupta et al., 2011a), which is 
used here. These combined schemes are calibrated 
with respect to the model ID convection equation and 



the contours of numerical amplification factor I G I and 
numerical group velocity V g N/C, are shown in Figs. 3(b) 
and 3(c), respectively. Despite providing superior 
discretization properties, one notices a range of high 
wavenumber components (kh > 2.40) displaying 
spurious dispersion, with V s n/C contours taking 
negative values. Such spurious upstream propagation 
gives rise to q-waves, which has been removed for a 
small Nc range by using a fifth order upwind filter 
(Sengupta et al, 2009), whose transfer function and 
effects on numerical properties are shown in Fig. 4. 

It is known that high accuracy methods also suffer 
from aliasing error for high Reynolds number 
computations. This has been controlled by using an 
adaptive multi-dimensional filter (Bhumkar and 
Sengupta, 2011), whose transfer function is shown in 
Fig. 5 for different filter coefficient values. Resultant 
method is calibrated by comparing experimental 
results (Fujino et al, 2003) and numerically-obtained 
time-averaged lift coefficient with angle of attack, for 
Re = 2.8 million in Fig. 6(a). Time histories of lift and 
drag coefficients in Fig. 6(b) and 6(c), show good 
match between computational and experimental time 
averages for Re = 10.3 million. Variation of RMS 
component of azimuthal disturbance component of 
velocity is shown in Fig. 6(d), showing flow transition 
to occur after 0.60c on the top surface and after 0.70c 
on the bottom surface. 

Numerical results have also been obtained for Re = 
10.3 million, at zero angle of attack, using a very fine 
grid with 5269 X 577 points, which are shown in Figs. 
7 to 9. Results show that due to high adverse pressure 
gradient, initially disturbances originate near the 
trailing edge on both top and bottom surfaces and 
with time these propagate upstream. Upstream 
movement of these disturbances is restricted up to x = 
0.55c on top surface and x = 0.70c on bottom surface of 
aerofoil. The top surface data matches with the 
experimental data in (Fujino et al, 2003), while there 
are no data available for transition on the bottom 
surface. 

Present exercise needs to be carried out with the 
application of realistic explicit background 
disturbances in wind tunnels to highlight the role of 
freestream turbulence in fixing the loads on aerofoil, 
specifically the drag coefficient. This would require 
not only the use of fine grids, but also proper 
description of free stream turbulence in the tunnel and 
its appropriate model. 
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